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FOREWORD 



For more than a decade, scientists from many different fields and countries have 
come together for the annual Optimization Days in Montreal to share their 
newest research results. This event is, in turn, hosted by various Departments of 
the four Universities in Montreal, viz. Concordia University, McGill University, 
University de Montreal and University du Quybec k Montryal. In 1986, the 
Optimization Days were hosted by the Ecole des Hautes Etudes Commerciales, the 
Business School affiliated with University de Montryal. In that year, for the first 
time, it was attempted to produce a collection of papers presented at that 
conference. The result is this volume. Its contents are divided into the following 
six sections: 

Mathematics of Operations Research and Global Optimization 

Linear and Combinatorial Programming 

Tours, Locations and Scheduling 

Dynamic Programming and Game Theory 

Control Theory 

Economic Models 

Each of the contributions in this volume has been refereed. This task was made 
easy by the many helping hands in the form of colleagues from Canada, Italy, 
Poland, U.S.A., Britain and The Netherlands. We wish to express our sincere 
thanks to them; without their help this volume could not have been completed. 

Additional thanks are due to our students Linda Mimeault and Ernest Leung and 
to Mrs. Halina Monkiewicz who typed the entire manuscript. 



Montreal, Canada 
October 1987 



H. A. Eiselt and G. Pederzoli 
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MATRIX DIFFERENTIAL EQUATIONS AND LYAPUNOV 
TRANSFORMATIONS 



John Jones Jr. 

Air Force Institute of Technology 
Dayton, OH, USA 



ABSTRACT 

The main purpose of this paper is to establish necessary conditions and 
sufficient conditions for the existence of a solution Q(t) of the Lyapunov matrix 
differential equation of the form 

= A(t) Q (t ) _ Q(t)B(t) + C(t), t 6 [7, «] 

where Q(t), , Q * (t), * are hounded for t e [7, S[. Such 

solutions are used to generate Lyapunov transformations 'which map one 
differential system to another system which has the same stability properties. 
Solutions of higher order Riccati matrix differential equations are used to uncouple 
stiff linear time- varying differential systems. 

1. INTRODUCTION 

Capital letters will denote n by n matrices whose elements are continuous 
functions of a real variable t for t e [7, ti\. C nXn (t) will denote the vector 
space of n by n matrices having elements which are continuous functions of t for 
t e [7, tf]. 

In model simulation theory, optimal control theory, and stability theory of 
systems of linear differential equations it is useful to be able to transform given 
systems to another whose stability properties are known or are more easily 
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obtainable. The basic problem to be considered here is to determine how to 
obtain such a transformation. 



Definition . A non-singular matrix Q — Q(t) is a Lyapunov transformation if it is 
continuously differentiable, Q(t), Q *(t), an< * ' 



are each bounded for 



t 6 [7, 6 ] where -00 <7<t<tf<oo. Given two differentiable systems 
— A(t)x and = B(t)^(t). They are said to be L-equivalent to each other 
if both have solutions for t € [ 7 , 0] and there exists a Lyapunov transformation 
Q = Q(t) such that the following pair of equations 



B(t) = Q V)A(t)Q(t) - Q V) ^ 
A(t) = P _1 (t)B(t)P(t) - P~V) ^ 



( 1 . 1 ) 



have solutions where P(t) = Q *(t) holds. Necessary conditions and sufficient 
conditions for the existence of a solution will be established. 



Sufficient conditions for which the stiff linear system of differential equations 
of the form 



dx(tl _ 

dt “ 



-B(t) 1 pm 

-C(t) I A(t) 

y 



x(t), 



t e [% tf] 



( 1 . 2 ) 



may be uncoupled will be established. Such uncoupling techniques are shown to 
be related to nonlinear Riccati matrix differential equations. 



Results obtained in this paper extend classical Floquet theory for linear 
time-varying differential systems having periodic matrices as a state matrix to 
more general time- varying linear differential systems. 

The first equation in (1.1) above implies that 

^P = A(t)Q(t) - Q(t)B(t) (1.3) 



and the second equation in (1.1) above implies that 




3 



^ - B(t)P(t) - P(t)A(t) 



(1.4) 



where P(t) = Q*(t)> 4 e [ 7 , $]. Upon multiplying on the left-hand side of 
equation (1.4) termwise by Q(t) and termwise on the right-hand side by Q(t) we 
get the following: 



Q(t) 



<m± 



dt 



-1 



Q(t) = Q(t)B(t) - A(t)Q(t) 

,-l 



(1.5) 



Now Q (t)Q(t) = I, 



d m 



dt 



+ dQM = 0> Mti = _ Q(t) 



-.-l 



Q(t) 



dQ(t) 

dt 



dt dt 

= Q(t)B(t) - A(t)Q(t), t e [ 7 , S\ 



Q -1 (t) ^ - 0, Q(t) ^ 

Q(t) . Thus using (1.5) it follows that 

( 1 . 6 ) 



=A(t)Q(t) - Q(t)B(t), t e[7, S\ 



(1.7) 



which is the same as equation (1.3). Therefore if equation (1.4) holds then so 
does equation (1.3). Similarly if equation (1.3) holds then equation (1.4) holds for 
t e [ 7 , tf] and x(t) = [Q(t)] y(t) maps = A(t)x(t) to - B(t)y(t). 



2. UNCOUPLING OF THE STIFF SYSTEM (1.2) 

In this section the problem of uncoupling stiff-linear differential systems of 
the form (1.2) will be considered. 



Theorem 1 Let the stiff-linear differential system (1.2) be given where A(t), B(t), 
D(t) e C nXn (t), t € [7, S\. If Q(t) is a bounded solution of the following Riccati 
matrix differential equation 

- A(t)Q(t) + Q(t)B(t) + C(t) + Q(t)D(t)Q(t), t e [ 7 , S\ (2.1) 

ax 

then there exists a Lyapunov transformation given by 
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( 2 . 2 ) 



( 2 . 3 ) 



y(t) 




x(t) = p 1 (t)x(t) 



such that the system (1.2) maps to the following system 




and conversely. 



t e [7, tf] 



( 2 . 4 ) 



( 2 . 5 ) 



Proof - Let Q(t) be a bounded solution of (2.1), then x(t) = P(t)y(t) where y(t) 
is given by (2.2). Upon differentiating x(t) = P(t)y(t) we obtain the following: 




Thus from (2.6) we see that 




and Anally 
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dy(t) _ 
dt ” 



-Q(t)B(t)-C(t) 


Q(t)D(t)+A(t)l 




' 0 


I 1 




0 


dQ(t)\ 

dt 


-B(t) 


D(t) 




I 


-Q(t) 


+ 


0 


0 




J 






) 






i 



Q(t)D(t) A(t) 


-Q(t)B(t)-C(t)-Q(t)D(t)Q(t)-A(t)Q(t) + ^ 


D(t) 


-B(t)-D(t)Q(t) 

; 



> y(t) 



y(t) 

( 2 . 8 ) 



Therefore if Q(t) satisfies (2.1) then (2.5) follows. 



3. EXAMPLE (THEOREM 1) 

Given the following continuous linear differential system: 



x(t) 



, 2 +t-l 








\ 


t+1 


1 


x(t) = 


-B(t) 


D(t) 


-1 


2 

t +3t+l 


-C(t) 


A(t) 

7 




t+1 > 





x(t) - A(t)x(t) 



(3.1) 



We need to solve the matrix Riccati differential equation (2.1) of the form 
2 2 

Q(t) - Q(t) + Q(t) + 1 + Q(t)iQ(t) (3.2) 

Thus 

Q(t) = Q 2 (t) + 2Q(t) + 1 = (Q(t) + l) 2 (3.3) 



has a constant solution Q(t) = -1, and so 



x(t) = P(t)y(t) = 



0 

1 



-Q(t) 



y(t) 



( . \ 



0 


1 


1 


1 

1 



\ ' / 



y(t) 



(3.4) 



and 



y(t) = p 1 (t)x(t) = 



f-i i\ 



V 1 0 / 



(3.5) 



Also 
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and x(t) = A(t)x(t) maps to y(t) — B(t)y(t) using the transformation (3.4). Also 
the characteristic roots of A(t) are obtained in (3.8). The form of 
y '(t) = B(t)y(t) given by (2.5) is given by (3.7) below: 




4. NECESSARY CONDITIONS FOR THE EXISTENCE OF A SOLUTION OF 
THE LYAPUNOV DIFFERENTIAL EQUATION 

Theorem 2 Let Q(t) be any non-singular solution of 

“iP = A(t)Q(t) - Q(t)B(t) + C(t) (4.1) 

where A(t), B(t), C(t) e C nXn (t) for t e [7, fl. Define R(t) by 




Then the following pair of matrices 




are similar for t e [7, £]. 

Proof - Let Q(t) be any non-singular solution of (4.1) above, then 




(4.4) 
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r -Qft)B(t)+C(t) 


Q(t)Q 1 (t)-A(t) 1 




a 


■ 1 


-B(t) 


q'VjqWQ -1 ^ 




i" 


-Q(t) 

/ 



r Q(t)Q V)-A(t) 


-Qft)Bft)+Cft)-Q(t) Aft)Q(t) 1 


^Q -1 (t)Q(t)Q -1 (t) 1 


-B(t)-Q _1 (t)Q(t) j 



and the matrices of (4.2) are similar. 



If Q(t) is any non-singular solution of (4.1) above then the pair of matrices 
in (4.3) are similar for t e [ 7 , £] and |R(t) — XI | is reducible to the product of 
at least one pair of polynomials f^(t, X) and g^(t, X), namely, the characteristic 

polynomials of the matrices Q(t)Q _1 (t)-A(t) and -B(t)-Q -1 (t)Q(t) respectively, 
having coefficients belonging to the ring of continuous functions of t e [ 7 , £] and 
such that f (t, X) is a degree a < n in X and g ft (t, X) is of degree a < n in 
X. Also tjt, R)g^(t, R) - 0. 

Now f (t, X)g^(t, X) is not necessarily the minimum polynomial satisfied by 
R(t) but is a divisor of |R(t) - XI| such that 



f a (t, Q(t)Q V)-A(t)) = 0; g/? (t, -Q V)Q(t)-B(t)) = 0 



(4.5) 



Definition . Polynomials f^(t, X) and g^(t, X) of degree a < n and of degree 
j} < n in X having continuous coefficients in t and such that f^(t, M * s 

a divisor of |R(t) - XI| and a multiple of the minimum polynomial satisfied by 
R(t) will be called an admissible class of polynomials and designated by OC(X) 
whenever the above properties of f^(t, X), g^(t, X) hold. 



Theorem 3 . Let f (t, X) e OC(X) such that tjt, Q(t)Q V) - A(t)) = 0 where 
Q(t) is any non-singular solution of equation (4.1) and 



f a (t, R(t)) = 



U M 
V N 



t e [ 7 , $] 



(4.6) 



and U, V, M, N are polynomials in the matrices A(t), B(t), C(t), 

-B(t)-Q -1 (t)Q(t), Q(t)Q _1 (t)-A(t), then Q(t) is a common solution of the pair of 
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equations 



Q(t)M + N =* 0 ; Q(t)U + V * 0 



(4.7) 



which may be written as follows: 

(Q(t) I) f a (t, R) = (Q(t) I) ^ - (0 0) (4.8) 

Proof : Let Q(t) be any non-singular solution of (4.1) then we have the following: 




and thus we see that equations (4.7) hold. 

Theorem 4 . Let g^(t, X) e OC(X) such that g^(t, -Q^ ^(t)Q^(t) + B(t)) = 
where Q^(t) is any non-singular solution of equation, and 



R) 



U mI , , „ 

Z — Z » 1 e [7, s\ 
V N, 

l / 



(4.10) 



where U, V, M, N are polynomials in the matrices A, B, C, -B(t)~Q *(t)Q(t) 
and -Q^(t)Qj(t) + B(t). Then Q^(t) is a common solution of the following 
pairs of equations: 



U - MQ 1 (t) ; V - N Q x (t) , t e [7, fl 



(4.11) 
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which may be written as follows: 



g^(t, R) 



( I ) 
QjW 



/- 

u 



— \ 
M 



N 



( I \ 
-Qj(t) 

J 



, t e [7, tf] 



(4.12) 



Proof. Similar as in the proof of theorem 3 above. 



5. SUFFICIENT CONDITIONS FOR THE EXISTENCE OF A SOLUTION OF 
EQUATION (4.1) 



Let 

R(t) 



-B(t) 


Q V)Q(t)Q W 


; yt, R) - 


U 


\ 

M 


C(t) 


-A(t) 


V 


N 




) 






/ 



(5.1) 



where f (t, X), g Jt, X) e OC(X) such that f (t, Q(t)Q *(t)-A(t)) = 0. 

OL fj OL 

g (t,-Q ^(t)Q(t)~B(t)) = 0 and Q(t) is any non-singular solution of equation 
(4.1). Now R(t) and f^(t, R) are commutative so that the following identities 

hold: 



(MC - UB = -BU + Q 1 QQ V ; MA - BM = UQ 1 QQ *Q *QQ *N 

-1 • -1 

|AV - VB = CU - NC ; CM - AN NA + VQ QQ 

for t e [7, $]. 



Theorem 5 . Let Q^(t) be any common non-singular solution of the pair of 
equations 



U - MQj(t) - 0 ; V - NQ t (t) - 0 (5.3) 

where N * exists for t e [7, 0]. Then Q^(t) is also a solution of (4.1) for 
t 6 [7, £]. 

Proof. Consider the following under the assumption that (5.1) holds and N * 
exists, then we have the following upon making use of (5.2) and (5.3). 
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0 = -A(V - NQj) (5.4) 

= (-AV) + (AN)Qj 

= (-CU + NC - VB) + (CM + NA - VQ" 1 QQj^Qj 
= (-CMQj + NC - VB) + CMQj + NAQj - VQ"^ 

= (NC - VB) + NAQj - VQ” 1 Qj 

- NC - NQjB + NAQ X - NQ^" 1 
= N(C - QjB + AQj - Qj) 

Now N * exists by assumption above and hence Q^(t) is also a non-singular 
solution of (4.1) for t e [7, 6\. 

Theorem 6 . Let Q(t) be any non-singular solution of the pair of equations given 
by (4.7), namely, QM + N = 0, QU + V = 0 where U * exists for t e [7, 6\. 
Then Q(t) is also a solution of (4.1) for t e [7, 5]. 

Proof . Consider the following under the assumption that (4.7) holds, then we 
have the following if U * exists by making use of (5.3) 

0 (QM + N)C 

QUB - QMC + QUB - NC (5.5) 

- Q(UB - MC) - QUB - NC 

= Q(BU - Q -1 QQ _1 V) - QUB - NC 

AQU + QBU + AQU - QUB - NC - QQ _1 V 

AQU + QBU - AV - NC + VB - QQ -1 V 
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- QBU - AQU - CU + QU 
= (QB - AQ - C + Q)U 

Now since U * exists by assumption above Q(t) is a non-singular solution of (4.1) 
for t e [7, S\. 

Theorem 7 . Let Q^(t) be any non-singular solution of the pair of equations given 
by equations below, namely, 

U = MQ X ; V = NQ X (5.6) 

where M * exists for t € [7, 6]. Then Q^(t) is also a solution of (4.1) for t e 

H *]• 

Proof . Consider the following by making use of (4.11) and (5.2) we have 

0 = (BU + MC - UB - Q^QQjQ^V) - (-MA + BM - Q" 1 Qj N + 
UQ" 1 Q Q^Qj 

= (BMQj + MC - MQjB-Q” 1 QjNQ) + (MAQj - BMQj + Q^QjQ^NQj 
MQj) 

= -M(-C + QjB - AQj + Qj) 

Now M * exists by hypotheses and hence Q^(t) is a non-singular solution of 
equation (4.12) for t e [7, £]. 

Theorem 8 . Let Q(t) be any non-sinular solution of the pair of equations QM + 
N = 0, QU + V = 0 where M * exists. Then Q(t) is also a non-singular 
solution of equation (4.1) for t e [7, 6\. 
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Proof. Consider the following by making nse of (4.7) and (5.2) then we have 
0 (QM + N)A 

Q(MA) - NA (5.8) 

Q(BM + UQ _1 Q Q -1 - Q~ 1 Q Q _1 N) - NA 

QBM - QUQ -1 Q Q -1 + Q Q _1 N - NA 

QBM - AQM + AQM + Q(-M) - NA - QUQ -1 Q Q -1 

= -QBM - AQM + AQM - QM - NA + VQ -1 Q Q _1 

QBM - AQM + AQM - QM + CM - AN 

QBM - AQM + AQM - QM + CM + AQM 

QBM + AQM - QM + CM 

= (-QB + AQ + C - Q)M 

Now M exists by hypotheses above and hence Q is a non-singular solution of 
(4.1) for t e [7, $]. 

6. EXAMPLE OF A LYAPUNOV MATRIX DIFFERENTIAL EQUATION 
Given the equation of the form 

A(t)Q(t) - Q(t)B(t) = ^ , t > 0 (6.1) 

where 




( 8 . 2 ) 
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Now 




where J[A(t)], J[B(t)] are the Jordan Canonical forms of A(t), B(t) respectively. 
Thus let 

c" 1 A(t)C 1 = J[A(t)] ; c" 1 B(t)C 2 = J[B(t)) (6.5) 

Using equations (6.1) and (6.5) to get 

[C^AMCj] [C- 1 Q(t)C 2 ] - [C^QWC^ [C~\t)G 2 ] - £ [C^Q^l (6.6) 
which can be solved for C^Q^cJ to finally find the solution of (6.2) to be 




(6.7) 
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It is also seen that , Q *(t), — - [Q *(t)] exist for t > 0 and thus 

Q(t) provides a Lyapunov transformation which maps = A(t)x(t) to = 

B(t)y(t) and both differential systems have the same stability properties for t > 
0. In this case both systems become unstable as t — ► +oo since the matrices 
A(t), B(t) become unbounded as t — ► +oo. 

SUMMARY 

The methods of approach used to establish necessary conditions and sufficient 
conditions for the existence of solutions to the Lyapunov matrix differential 
equation can be extended to nonlinear matrix Riccati type differential equations. 
Such results will appear elsewhere. For a detailed early treatment of Lyapunov’s 
matrix differential equation see references A. M. Lyapunov (1), J. F. P. Martin 
(2), V. V. Nemystkii and V. V. Stepanov (3) below. Results obtained in this 
work are useful in treating linear differential systems containing parameters in 
addition to uncoupling stiff linear differential systems and determining stability 
properties of linear differential systems. 
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THEORY AND METHODS FOR GLOBAL OPTIMIZATION 
INTEGRAL APPROACH 

Quan Zheng* 

Shanghai Institute of Applied Mathematics and Computation 
Department of Mathematics 
Shanghai University of Science and Technology 
Shanghai, People’s Republic of China 



1. INTRODUCTION 

1.1 Statement of the Problems 

Let X be a Hausdorff topological space, S c X a closed 
f = X — ► R a real-valued function. The problem considered here is to 
infimum of f over S, 

c = inf f(x) 
xes 

and the set of all global minima: 

H = {x | f(x) = c, x e s} 

We assume in this paper: 



— AN 



set and 
find the 

(i.i) 



( 1 . 2 ) 



This research was supported by the Science Foundation of the Academy of 
Sciences of China. The paper was finalized in May 1986 while the author was 
visiting professor at the University of Quebec in Montreal, Canada. 
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(Aj) f is continuous: 

(A ) There is a e R' such that the level set 

u 

H a - ( x I f M ^ “> ( L3 ) 

is compact and n s ^ 0. 

Thus the problem (1.1) becomes to find 

c = min f(x) = min f(x) (1.4) 

xes xeH nS 

a 

and the set of all global minima H is non empty. 

Most of the current theory and methods for dealing with this minimization 
problem are gradient-based. They are difficult to use for the solution of a global 
optimization problem. In this paper the integral-based theory and methods are 
presented. The advantages of this approach are that it can be used to solve 
global optimization problem with nonconvex, nonsmooth and even discontinuous 
functions. [1] 

We first introduce the concept of Q-measure defined on a normal topological 
space which has the property that each non empty open set has its positive 

measure. With the help of this tool the mean value of a function over its level 
set is introduced and investigated. Necessary and sufficient conditions for global 
optimality, in terms of the behavior of the above mean value, are then derived. 
We then extend the mean value definition to the case of variance as well as 

higher moments of a function and obtain the corresponding characterizations of 
global optimality for these cases. The constrained cases are examined by the 

rejection and the reduction methods. A modification of the penalty method for 

the mean value as well as for the variance and higher moment characterizations 
are discussed. The relationship between our approach and the standard, gradient 
based Kuhn- Tucker theory, convex analysis and Clarke’s nonsmooth cases is 
explored. With a suitable definition of continuity with discreet topology, integer 
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and mixed programming problems are considered. 

We then describe the mean value-level set method and prove its convergence 
(which does not depend on the choice of initial data) and the limited influence of 
the error propagation. The corresponding rejection and reduction algorithms, and 
the global version of the nonsequential penalty algorithm are discussed for 
constrained minimization problems. 

The optimality conditions and algorithms for the integral global optimization 
can be implemented by the use of properly designed Monte Carlo technique. 
Numerical tests and applications show that this approach is effective, in 
particular, for engineering problems. 

1.2 Q-measure spaces 

Let X be a normal topological space, A a tr-field of subsets of X and a 
measure on A. 

Definition 1.1 A measure space (X, fl, p) is said to be Q-measure space if 

(1) fl is a Borel field; 

(2) /i(G) > 0 for each non empty open set G e fl; 

(3) n(K) < +oo for each compact set K 6 0 

The Lebesque measure space in R n , (R n , /?, /i) is a Q-measure space. The 
non-degenerate Gaussian measure on a separable Hilbert space H, (H, fl, /i) is 
also a Q-measure space. Let X = {x^ ..., x^, ...} with discreet topology and 
X 

fl = 2^ . For each set A e fl we define j*(A) = {][]ai|x. e A}, where ai > 0, 

i 

00 

i = 1, 2, ..., and a i = 1. Then (X, fl, /i) is a Q-measure space too. A 

i=l 

specific measure space can be applied to solve a specific minimization problem. 
In the following discussion we always suppose that we have a Q-measure space. 
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Before dealing with general optimality condition of global minimization, we 
first give a lemma which can be regarded as a sufficient global optimality 
condition. 

Definition 1.2 A subset F of a topological space X is said to be robust if 

cl (int F) = clF (1.5) 

Lemma 1.1 . Let (X, Q, /i) be a Q-measure space and S c X a closed robust 
set. Suppose that the intersection of the level set = {x | f(x) < c} and S is 
non empty. If the measure 

#*(H. nS) = 0, (1.6) 

Then c is the global minimum value and n S is the set of all global 
minimizers of f over S. 

Corollary , if c > c = min f(x), then 

xeS 

fi { He n S) > 0 (1.7) 

This corollary is useful to define the concept of the mean value. 

2. INTEGRAL CHARACTERIZATION OF GLOBAL OPTIMALITY 
[2, 3] 

2.1 Mean Value Conditions [4] 

We begin with the unconstrained problem. Let c = inf f(x) 

xeX 

Definition 2.1 Suppose c > c, we define 

M(f ’ c) = (ife ^ He f(x) d/i 

to be the mean value of f over its level set = { x | f(x) < c} 



( 2 . 1 ) 
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It is easy to verify the following statements: 



(1) Homogeneity M (af, ac) = aM(f, c), for constant a > 0 and c > c; 

(2) Translation M(f+a, c+a) = M(f, c) + a, for constant a and c > c; 



(3) For c > c, M(f, c) < c; 

(4) Monotonicity : If c^ > c^ > c, then M(f, c^) M(f, c^); 

(5) If {C^} is a decreasing sequence whose limit is c > c, then 

M(f, c) — lim M(f, c^) 

C, — *c 
k 

When c — c, /i(H_) might equal zero. We then extend the definition of 
mean value by a limit process. 



Definition 2.2 Let c > c and {c^} be a decreasing sequence whose limit is c, 
then the mean value M(f, c) is defined to be 



M(f, c) - 



hm /i(Hc ) f He dfi 

c. —vc R' R 



( 2 . 2 ) 



We can prove that (2.2) is well defined and independent of the choice of 
{C. }. By the statement (5), definition (2.2) is consistent with (2.1). All the 

Jv 

statements (1) — (5) hold for this limit- based definition. 



Theorem 2.1 (Mean Value Conditions). For the problem (1.1) under assumptions 
(A ) and (AJ, a point x is a global minimizer with i = f(x) as the corresponding 

X M 

global minimum value of f if and only if 



M(f, c) > c, for all c > c 



(2.3) 



or 



M(f, c) = c. 



(2.4) 
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Example . f(x) — |x| a with a > 0. H c = x{x| |x| a < 0} * [-c^ a , c^ a ] for 
c > 0. M(f, c) = 1/1+a c. Applying (2.4), we have 1/1+ac = c which implies 
c = 0 and H _ — {o}. 

2.2 Variance and Higher Moment Conditions [5] 

We now go a step further to introduce the concepts of variance and higher 
moments. 



Definition 2.3 Suppose c > c. We define the variance 

T(f ’ c) = dn /h {f(x) ‘ M(f> c))2<1/i ( 2 - 5 ) 

^ c 7 c 
and the modified variance 

V l (f) c) - MHM /h (f(x) " C)2d/i 
v c 7 c 

of f over its level set H . 

c 

The variance of a function has the following properties: 

(1) Positivity V(f, c) > 0; 

2 

(2) Second-degree homogeneity: V(af, ac) = a V(f, c), for a > 0 and c > c; 

(3) V(f+a, c+a) = V(f, c) for c > c; 

(4) Vjff, c) - V(f, c) + (M(f, c) - c) 2 ; 

(5) If{c fc > is a decreasing sequence which tends to c > c, 
then V(f, c) - lim V(f, c k ) 

c i 

k— >c 

Definitions (2.5) and (2.6) can be extended for c > c by a limit process. 
And all of the statements (l)-(5) hold for these limit-based definitions. 
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Theorem 2.2 (The variance conditions). For the problem (1.1) under assumptions 
(Aj) and (A^), a point x is a global minimizer with c = f(x) as the 
corresponding global minimum value if and only if 



V(f, c) - 0 (2.7) 

or 

Vj(f, C) = 0 (2.8) 

Example The criteria for Nelder-Mead Simplex Method of nonlinear minimization: 



■ E ( f ( x i) - f ( x c )) 2 < £ . x c = ( x i + ••• + x n )/ n 

i=l 



1 n 2 

- Y (f(x.) - f ) < e, f = (f(x ), ..., f(x )) 

n v x v max' max u r v n" 



(2.9) 

( 2 . 10 ) 



can be regarded as an approximation forms for the variance conditions, if f is 
convex. 



Definition 2.4 Suppose that m is a positive integer, a is a constant and c > c. 
We define 

Mm(f, c; a) = (f(x) - a) m dyi (2.11) 

to be the m-th moment of f over its level set H centered at a. 

c 

By a limit process the definition can be extended also for c > c . 

Theorem 2.3 (The Odd Moment Conditions). For problem (1.1) under 
assumptions (Al) and (A2), a point x is a global minimizer with c = f(x) as the 
corresponding global minimum value if and only if for some positive integer m we 
have 

M 2m-1 ( f ’ c; °) - for c > c (2.12) 

or 
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M 2m-1 (f ' °) “ ^ 2m 1 ( 2 - 13 ) 

Theorem 2.4 (The Even Moment Conditions). With respect to the above 
problem, a point x is a global minimizer with c = f(x) as the corresponding 

global minimum value if and only if for some positive integer m we have 

M 2m (f, c; M(f, c)) = 0 (2.14) 

Theorem 2.5 (The Hogher Movement Conditions). A point x is a global 
minimizer for the above Theorem, and c = f(x) is the corresponding global 
minimum value if and only if for some positive integer m we have 

M m (f, c ; c) = 0 (2.15) 

2.3 Rejection and Reduction Conditions [6, 7] 

The integral chracterizations of global optimality have a unified form both 

for unconstrained and constrained problems. Suppose that the admissible set S is 

robust, then we can construct a Q-measure such that = A*(A n s), and 

(X n s, fls, /*s) is called the rejection measure space, where fls = {A n s/ A e 

0}. We can also construct the Q-reduction measure space on a manifold of X, 
(x) 

say, S = {x hr } = 0, i = 1, ..., r}. We now easily define rejection and 
reduction mean value, variance and higher moments M(f, c; s), V(f, c; S), Vl(f, c; 
S) and Mm(f c; a; S), respectively. For constrained problems we also have 

similar characterizations of global optimality. 

Theorem 2.5 With respect to the constrained minimization problem 

min f(x) 
xes 

Under assumptions (Al) and (A2), the following statements are equivalent: 
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(1) x g S is a global minimiser and c = f(x) is the corresponding global 
minimum value; 

(2) M(f, i; S) > c, for c > c; 

(3) M(f, i; S) - c; 

(4) V(f, i; S) - 0; 

(5) Vl(f, i; S) = 0 

2m 1 

(6) M ^ (f, c; 0; S) > (c) , for c > c and some positive integer m; 

J 

(7) ^(f, c; 0; S) = (c) , for some positive integer m; 

(8) M^^f, ^ S) = 0, for some positive integer m; 

(9) M^(f, c; c; S) = 0, for some positive integer m. 

2.4 Penalty Global Optimality Conditions [8] 

Suppose that X is a metric space and that S is a robust closed set. A 
continuous function p on X is a penalty function for S if (i) p(x) > 0 for all 
x e X and (ii) p(x) = 0 and only if x e S. 

Suppose that {c^} is a decreasing sequence which tends to c > c as k — ► oo 
and {ak} is a positive increasing sequence which tends to infinity. Let 

H k = {x | f(x) + o k p(x) < c k }, k = 1, 2, ... (2.16) 

Under assumptions (Al) and (A2) we can prove 

lim H = H n S (2.17) 

i X c 

k — ►oo 

“ d k'^= \ n * )d " ~ MH^s) H « " S ^ 



(218) 
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The limits are independent of the choices of {c^} and {ak}. (2.17) 

and (2.18) can be also extended for c > c. Thus, the penalty mean value can 
be defined as follows. 

Definition 2.5 Suppose c > c. The limit 

M(f, c; p) = lim -rg-r / R f(x)d/i (2.19) 

k— >oo ^ k' k 

is called the penalty mean value of f over its level set with respect to the 
penalty function p. 

We can define the penalty variance and higher moments and prove the 
penalty global optimality conditions analogous to Theorem 2.5. It is quite nature 
to define penalty mean value as 

M'(f, c; p) = lim "Vg" : / H (f(x) + ak p(x)) djt (2.20) 

In fact, we can prove 

M'(f, c; p) = M(f, c; p) for c > c (2.21) 

The same can be done for the penalty variance and higher moments. 

2.5 Convex and Nonsmooth Cases [9] 

Let X be a locally convex topological space, S be a convex set of X with 
int(S) 0, and f be a convex function on X. Consider the convex minimization 

problem: 

min f(x) (2.22) 

xeS 

The following lemma is the analogue of the separation theorem in convex 
analysis. 
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Lemma 2.1 A point x € S is a global minimizer of f over S with c — f(x) a 8 
the corresponding global minimum value if and only if there exists a subgradient 
£ € d f(x) such that 

M(f^, c; S) > c, for all c > c (2.23) 

where 

f ^ = c + <£, x - x> (2.24) 

The condition (2.23) is equivalent to the statement that the linear functional 
<£, x - x> is a supporting functional of S at x. 

Let X be a Banach space and f be locally Lipschitzvan. The following 
lemma provides the relationship between minima of a function and of its 
generalized directional derivative. 

Lemma 2.2 If x is a local minimizer of f, then 0 is a global minimizer of the 
generalized directional derivative f°(x; x). 

Now the global optimality conditions can be applied to f°(x; x). Since f°(x; 
x) is convex, we then reduce the problem to a convex one, and the following 
proposition can be proved. 

Proposition Suppose that x e S is a local minimizer of f over S, then 0 is a 
global minimizer of f°(x; x) over the cone T(x, S) of tangents of S at x (or over 
any nonempty closed convex cone T^ included in T(x; S)). 

2.6 Integer and Mixed Programming [10] 

Let X = x^, ..., x^, ... and t be the discrete topology on X. Let S be a 
subset of X and f be a real-valued function. Then each function f on topological 
space (X, T) is continuous. 

X 00 

Let ft = 2 and ai > 0, i = 1, 2, ... with ai < oo. We define a 

i=l 



measure fjt on by 
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m 

where 


= Y a. V A e ft 

•i 1 

ieJ 


(2.25) 


J = 


{i | x. e A}. 


(2.26) 



Then {X, Q, p} is a Q-measure space. We label 



I = {i | x. e H c n S} and p = /i(H n S) — ai (2.27) 

iel 

Definition 2.6 Suppose that c > c = min f(x) (under assumption (A 9 )). We 

xgS 1 

define the mean value, variance and higher moments of f over its level set H n 

c 

S, respectively, by 



M(f, c; S) = y aif(x. / ft c; 
iel 


(2.28) 


V(f, c; S) = y ai[f(x.) - M(f, c; S)] 2 / /xc; 
iel 


(2.29) 


V (f, c; S) = £ ai[f(x ) - c] 2 / /ie; 
iel 


(2.30) 


Mm(f, c; u; S) = y ai[f(x.) - a] m / fic. 
iel 1 


(2.31) 


With the help of these definitions, the theorem 


analogous to Theorem 2.5 



holds for the integer programming. 

Two kinds of mixed programming problems are considered here: 

(1) The case of product spaces, X =» I x Z , where I is discrete and Z is a 
normal space. We define a product Borel field to be 



Q = {A x B | A g 0, and B g Q2} 



(2.32) 



and the measure is given by 
#*(A X B) = #* 4 (A) • /x 2 (B), 



(2.33) 
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where the two Q-measure space (I, Q^, p^) and (Z, p^) are described above. 

The new measure space is also a Q-measure space. 

(2) The case of union spaces, X = I u Z, and I n Z = 0. An open set in X 
is the union of an open subset of I and that of Z. The union Borel field is 
given by 

fl = {A u B | A € 0, and B € n 2 } (2.34) 

and the measure is given by 
/i(A u B) = Mj(A) + h 2 {A), 

where the two Q-measure space (I, Q^, p^) and (Z, fi^, p^j are described above. 
The new measure space is a Q-measure space too. 

It is not difficult to define mean value, variance and higher moments and 
deduce the characterizations of global optimality for these mixed programming 
problems. 

3. THEORETICAL ALGORITHMS AND THEIR IMPLEMENTATION 

3.1 M-L Algorithm (Mean Value — Level Set) [11-13] 

Due to the fact that both unconstrained and constrained problems have the 
same characterization of global optimality we will only describe M-L algorithm for 
the unconstrained problem (under assumptions (Al) and (A2)). 

M-L Algorithm 

Step 0 : Take and e > 0; k = 0; c q := f(x Q ) (> c = minf(x)); 

Hc o - f* I *w - V 

C k + 1 - “< f - "k); H ' k+ 1 I <w s W> 



Step 1 : 
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Step 2 : V := Vj(fj, c^); If V > e then k := k + 1 and go to Step 1 

otherwise go to Step 3; 

Step 3 : c «— H «— Hc^^; Stop. 

Let e = 0 in the above algorithm, then we get two monotonous sequences: 



c > c, > ... > c, > c, , > ... 

o 1 k - k+1 - 


(3.1) 


and 




He > He > ... > H , > ... > He, , > ... 

o 1 ck “ - k+1 “ 


(3.2) 



which are bounded below so the limits exist. Let 



c = lim c, and H == lim H , = H H 
k— >oo * k >oo CJC k=l Ck 

Taking the limit for c^^ = M(f, c^), we have 

c = lim (k + 1) = lim M(f, c k ) = M(f, c) 

k— >oo k— >oo 

By the mean value condition (2.4), we have proved 



(3.3) 



(3.4) 



Theorem 3.1 The limit c in (3.4) is the global minimum value of f and H in 
(3.3) is the set of all global minimizers of f. 

The convergence of the M-L algorithm does not depend on the initial choice 
of x^. If f(x) is not a constant we can prove, 

'k =■ Vl : "(’’ck) 2 '•< H c k + l>' ,0r k - *• *■ ■ ( 3 - 6 ) 

The important fact is that the errors in each step of the algorithm will not 
be accumualted. 

Suppose we calculate c^ = M(f, c^) with error A^, and get d^ — c^ + A , 
and then calculate c^ = M(f, d^) with the error A^ and get d^ = + A^, 

and so forth. Let d = lim d, . By the mean value condition, we can easily 

k— >oo 
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prove. 

Theorem 3.2 d is the global minimum value of f if and only if 

lim A — * 0 (3.6) 

k->oo 

3.2 Nonsequential Penalty Algorithm 

Penalty algorithms are useful to solve constrained problems when the 
rejection algorithm is not efficient. A nonsequential penalty algorithm, different 
from sequential unconstrained minimization technique, is proposed, see [14]. 

From 3.1, we realize that the method for finding global minimizers requires 

the computation of a sequence of mean values and a sequence of level sets. 

Finding a mean value is tantamount to computing an integral of a function with 

several variables. The determination of a level set is, in general, more involved. 

But accuracy at earlier steps is not generally required by Theorem 3.2. This 

suggests that a Monte Carlo based technique of finding global minimizers is 

appropriate. The error by a Monte Carlo method is proportional to a / y/t, 

2 

where t is the number of sampels and or is the variance of sample distribution. 
2 

Since <r tends to zero as the mean value goes to global minimum value 
(Variance condition), so the Monte Carlo approximation becomes more accurate 

near global minimum value even though the number t of samples may be not 
very large. 

3.3 Monte Carlo Implementation of a Simple Model 

We first consider a box constrained minimization problem in R^: 

min f(x), D = {x = (x\ ..., x 11 ) | < x < b*, i = 1, ..., n} (3.7) 

xeD 

with a unique global minimizer x. Let be the smallest cuboid which contains 

the set n D, for k = 1, 2, ... : 
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D fe - {x = (x\ x n ) | ^ < x 1 < b\ i = 1, n}, (3.8) 

where 

a k = inf ^ I x “ C* 1 ' ^ e H ck n D> 

= sup{x 1 I X = (x 1 , x n ) e Hc k n D}, 

Then we have 

c = min f(x) = min f(x) (3.9) 

xeD xeD. 

k 

and we can prove that 

H = {x} = n \ (310) 

k=l 

Instead of M(f, c^) and V^(f, c^) in M-L algorithm we take M(f, c^; D^) and 
V^f, c^.; D^) at each iteration. The Monte Carlo implementation of this simple 
model is described as follows: 

(1) Approximating Hco and M(f, co) Let £ = (£^ ..., £*) be an independent 
n-tdple random number which is uniformly distributed on [0.1]". Let 

x 1 = a* + (b* - a*) • £*, i = 1, ..., n (3.11) 

Then x = (x\ ..., x 11 ) is uniformly distributed on D q = D. 

Take Km samples: x^, x^, ..., x^ m and compute the values f(x^), j = 1, ..., 

km. Comparing the values of the function at each point, we obtain a set W of 
sample points which contains t points corresponding to the t smallest function 
values FV[j], j = 1, 2, ..., t ordered by their values, i.e., 

FV[1] > FV[2] > ... > FV[t] (3.12) 

The set W is called the acceptance set which can be regarded as an 
approximation to the level set with c^ =* FV[1], the largest value of {FV[j]}. 
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Clearly, f(x) < c q for x e W. Also, the mean value of f over can be 
approximated by the mean value of {FV[j]}. 

c - M(f, c q ) » [F[l] + ... + FV[t]]/t (3.13) 

(2) Generating a New Domain by W. The new cuboid domain of dimension n 

D1 - {x - (x 1 , ..., x n ) | a^ < x 1 < b^, i - 1, ..., n} (3.14) 

can be generated statistically. The following procedure is proposed. Suppose that 
the random samples in W are T^ ..., T^. 

Let 



tJ), c\ = max^, .... T^), i - 1, .... n, (3.15) 

where 

T. - (T 1 , T n ), j - 1, t. 

We use 

a 1 = a l o - (<Tj - ^ | (t - 1) and (l - o\ + {o\ - o q / (t - 1) (3.16) 

as estimators to generate a^ and b^, i = 1, ..., n; a and $ in (3.16) are 

unbiased estimators of the end points of an interval, if t uniformly distributed 
random samples are taken from the interval. 

(3) Continuing the Iterative Process The samples are now taken in the new 
domain D^. Consider a random point x, where 

X 1 = a^ + (bj - a'j) ■ i, i = 1, .... n 

compute f(x). If f(x) > FV[1] then drop it. Otherwise, reconstruct {FV[j]} and 
W such that the new {FV[j]} is made up of the t best function value obtained 

so far. The acceptance set W is modified accordingly. Repeating this procedure 

until FV[1] < C^, we obtain new FV and W. 
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Continuing this process gives us a decreasing sequence of mean values {C^) 
and a sequence of cuboids {D^} of dimension n. 

(4) Iterative Solution At each iteration k, the smallest value FV[t] in the set 
{FV[j]} and the corresponding point in W can be regarded as an iterative 
solution. 

(5) Convergence Criterion The modified variance VI of {FV[j]} which is given 
by 

1 * 2 

Vj = t -rr E ( FV tj] - Fvfin, 

j=2 

can be regarded as an approximation to Vl(f, c^). If VI is less than the given 
precision e, then the iterative process terminates, and the current iterative solution 
in (4) serves as an estimate of the global minimum value and of a global 
minimizer. 



Under suitable assumptions we can prove (see [13]) 



Theorem 3.4 The number of computations of a function f for capturing the 
global minimizer in a small cuboid of volume 6 from an initial cuboid of unit 
volume has the following asymptotic bound as ^ goes to zero. 



NL < C, in 
f ~ f 







(3.18) 



where C^ is a constant depending on f but independent of 6 . 

3.4 Remark on Other Models 

The technique in 3.3 can be extended to the case when the function f has 
multiple global minima. The search domain at the k-th iteration can be 
decomposed into unions of several cuboids of dimension n: 
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k k 

\ = U D-. 

j=l J 

so that each smaller cuboid D. can be treated individually as in 3.3. Usually we 
assume that r^ < m for each k with m as a positive integer given in advance. 

The rejection method is used to generate a probability on S n by the 
Monte Carlo technique. The reduction method is usually used in the equality 
constrainead problems by generating a probability on a manifold. It is not 
difficult to generate a random number £ with a given discrete distribution which 
can be used for integer and mixed programming problems. 

4. Numerical Tests and Applications 

Numerical tests have been done covering the area of unconstrained, equality 
constrained, and inequality constrained minimization problems. Here we list the 
tests for the unconstrained problem, the problem with multiple global minimizers, 
the nonlinear integer programming problem and a discontinuous minimization 
problem. 

(1) Unconstrained Minimization [16] 

minimize f(x) 
where 

f M = ^ {ksin(IIx 1 ) + Y, [ x j - A ) 2 (! + ksinfllx. + 1)) + (x q - A) 2 ] 

k=l 

k = 1 and A =* 10 

D ** {x = (x\ ... x 11 ) | - 10 < x 1 < 10, i = 1, ..., n} 

The following table gives the numbers of iterations Ni and the amounts of 
function computations N^ corresponding to n = 5, 10, 20, 50, when minimal 

function values are less than 10 ^ (with IBM PC): 




34 



n 


5 


10 


20 


50 


Ni 


61 


113 


231 


457 


N f 


2864 


6800 


16995 


53879 



(2) Multiple Global Minimizers The following table is taken from Zhou’s master 
thesis [15], where the algorithm was tested on the functions GOLDPR and RCOS: 

GOLDPR = f(x, y) = [1 + (x + y + 1)^ (19 - 14x + 3x^ - 14y + 6xy + 3y^)] 

[30 + (2x - 3y) 2 (18 - 32x + 12x 2 + 48y - 36xy + 27y 2 )] 

D 0 = {(x, y) I -2 < x < 2, -2 < y < 2} 

RCOS = f(x, y) = a(y - bx 2 + cx - d) 2 + e(l - f) cosx + e 

a = 1, b = 5.1 • (4II 2 ), c = 5/n, d = 6, e = 10, f = l/8n 
— {(x, y) | -5 < x < 10, 0 < y < 15} 

Zhou’s results are compared with those of the clustering method [16]: 
the table which goes back to original 

(3) Nonlinear Integer Programming Find the global minimizer of f over S, 
where f = GOLDPR and 

S - {(x, y) | x, y - O.Olj, j 200, -199, ..., 198, 199, 200} 

After twelve iterations, the variance VI is reduced from 0.77.10** to 0.0 and 
x = (0.0, —1.0), c = 3.0, Nc = 621. The function has four local minimizers if 
it is considered on 

D = {(x, y) | -2 < x < 2, -2 < y < 2} 
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(4) Discontinuous Minimization 

Find the minimum of f over D, where 



n + 



f(x) 



(E 

i=l 



I x.|) J + Sgn[8in(( £ 



I xjl) 5 ~ 0.5)] 



x ^ 0 
x = 0 



This is a discontinuous function with countable many discontinuous 
hypersurfaces as well as countably many local minima. The function f has an 
essential discontinuity at the origin which is also the unique global minimizer. 
With n = 5, after 100 iterations and N^ = 4962, we obtained on a personal 
computer the global minimizer. 



a = (-8.94.10 14 , 7.21 • 10 4 , -3.01 • 10 13 , 8.57 • 10 13 , 3.05 • 10 14 ) 



Integral global optimization algorithms have been successfully applied to many 
industrial and control problems such as automatic design of optical thin-film 
systems [17], automatic generation of prototype lenses [18], computer aided design 
of equalizers [19, 20], optimal design of microwave filters [21], optimal design of 
the turbine wheel, optimal design of speed reducers [22] and nonlinear of servation 
of dynamic systems [23]. 
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ABSTRACT 

A set-monotonie non-gradient algorithm is proposed for finding global minima 
of general non-convex mathematical programming problems. The algorithm is 
based on the Cubic Algorithm /l/ equipped with a semi-certain distinction 
operator and the marginal comparison constant generator. An improved version 
of the algorithm is presented as compared with /2 /. 

1. INTRODUCTION 

We consider the problem 

min f(x), x e X c R n (1.1) 



of finding the value 



p° = min f(x) 
xeX 

and the set 



( 1 . 2 ) 



X° = {x| f(x) = p°, x £ X}, (1.3) 

where X is a compact robust set which may be non-convex and non-connected. 
For example, X may consist of a collection of closed tori, balls and cubes that 
may intersect one another. 
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By definition, a set Y is called robust, if the closure of its interior coincides 
with its closure: cl int Y — cl Y. Note that Y may be open or closed, or 
neither open nor closed. A closed robust set X always has a non-empty interior. 
An open set X is always robust, its closure is denoted by X. 

Lemma . If X c R n is robust and compact, f: R n — ► R is continuous over 

X and c = const, then the set 

Y = {x | f(x) < c, x e X} (1.4) 

is either empty, or robust for every c. 

Introduce a circumscribed (non-strictly) closed cube C such that 

X c C c R n (1.5) 



and let c be the length of the edge of C. 



Hypothesis . The cost function f: R n — ► R is assumed to be a continuous 

single-valued computable procedure defined and Lipschitzian on C, that is 

|f(x) - f(x / ) | < L||x - x'||; L = const > 0; x, x' e C (1.6) 

Take a partition integer N > 2 and consider deletion constants: 



r 

m 



Lcy/n 

N m 



m = 1, 2, ... 



(1.7) 



The notation J|.|| stands for Euclidian norm. 



2. THE SEMI-CERTAIN DISTINCTION OPERATOR 

Definition 1 . Given X, C, a distinction operator is defined by the binary 
function: 



4> = $(X, C) 



0, if X n C = <f> for sure 



1, otherwise 



( 2 . 1 ) 
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Here the quantifier “for sure* has nothing to do with probability (not to be 
confused with “almost sure”), also the term “otherwise” does not always mean 
X n C - 4>. 

In some cases it is easy to construct distinction operators. 

Example 1 . Suppose that C c R n is a closed cube with the edge c > 0 
and X c R n is a finite union of closed balls: 

X = {x| ||x - a^|| 2 < b 2 , b. > 0, 1 < j < k} . (2.2) 

For this case a distinction operator is given by the inequality: 

II* - a jll > dj, z e C, j = 1, k (2.3) 

where 



d. > b. + c^n, (2.4) 

with the understanding that <£ = 0, if (2.3) holds for any z c C fixed in 
advance, and <£ =* 1 otherwise. Distinction operator is not unique and for the 
above case the one given by the equality in (2.4) may be considered as “better” 
one. 

Example 2 . Let X - [-2, 2], C = [2, 3], = [2.5, 3.5], C 3 = [3, 4] 

and = [4, 5]. Then we have = $ 3 = 1 and = 0 meaning 

that only X n C. = </> for sure. In fact, X n CL ^ X n CL — X n CL 
4 1 

= <f> but all three not “for sure”. Indeed, this example is a special case of (2.2), 
(2.3), (2.4) with d > 3, n = 1, k = 1, a^ = 0, b^ = 2, c = 1 for all three 
cubes. Our distinction operator is, thus, as follows: 




1, otherwise 



(2.5) 

( 2 . 6 ) 
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Let d = 3 and take z^ =* 3 e C^. For this z^ we have (2.6): 0^ * 1 

since (2.5) is not satisfied. If we took d < 3, then (2.5) would be satisfied 
yielding 0^ = 0, in contradiction with the fact that X n «= {2} <j>. This 

illustrates that inequalities (2.3), (2.4) are unimprovable . 

Now, take z^ = 3 e yielding as before $ = 1 despite the fact that, 

had we taken z^ = 3.5 e C^, we would have (2.5) satisfied for this particular 

choice. The point here is that the statement “X n C. = ^ for sure” means that 
the empty intersection is to be established by the check based on one single 

point in C. arbitrarily fixed in advance which is essential for the operation of the 
Beta- Algorithm. Naturally, such an operator cannot be fully certain and this is 

the sense of the word “otherwise” in (2.1) comprising all cases in which the 

non-intersection cannot be established by the check of a single point z. e C.. In 

the example, those are cases for CL, C , C each of which has at least one 

1 A O 

point, e.g., z = 3 for all three cubes, violating (2.5). 

The situation is different for CL. Here, whatever z € CL, we always have 

4 4 

|z| > 3, thus, 0^ = 0. 

Let us now take d = 4 > 3 in (2.5). In this case 0^ = 1 (for z *= 4 e 

the inequality |z| > d is violated), however, for = [4.5, 5.5] we have 

O c 0, meaning that with such d the operator is still working, but it is poorer, 
having a greater uncertainty band. 

From the above examples we can see that the distinction operator generates 
a compact set Y containing X and such that, given any fixed z e C, we have 
X n C = 0, for sure, if z ^ Y, denoted in (2.1) as 0 = 0; otherwise, 0 = 1 
that corresponds to two possibilities: z e X meaning X n C ^ and z e Y - 

X in which case nothing is known about the intersection. For a “good” 
distinction operator the sets Y and X are congruent and such that for the cubic 
set C c R n with the edge c > 0 the boundaries dY and <9X are c>/n- equidistant 
surfaces, yielding the uniform uncertainty band V = Y — X. 
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Let us consider variable sets X, C parametrized by a parameter m c [0, oo): 
X m = X(m), C m = C(m) (2.7) 

where X(0) = X, C(0) = C and all C m are closed cubes with edges of the 
length 

c m = c/N m , m = 0, 1, 2, ... (2.8) 

In this case the uncertainty band D will depend on m : = P(m) = - 

X . 

m 

Suppose X - c X (m = 0, 1, ...) and that the intersection of all X is 
m+1 m v ’ ' m 

non-empty, which allows us to define the limit: 

oo 

lim X = f| X = X° ^ 0, clearly X° c X. (2.9) 

m — >00 m=0 

Suppose also that the limit lim Y is defined. 

m— *oo 

Definition 2 . A semi-certain distinction operator <£(m) ^(X^, ® m ) * s 

called precise, if 

lim D - <t>, (2.10) 

m— *oo 

or, equivalently: 

lim Y = lim X = X° . (2.11) 

mm 
m — >00 m— ►oo 

It means, of course, that the boundaries dY and dX tend to coincide as m —► 

mm 

oo. Clearly, the parameter m does not have to be discrete. One may choose a 
continuous parameter, for example, the length c of the edge of a variable cube C. 
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Precise distinction operators are not unique. For instance, suppose that 
— X = const given by (2.2). If we take in (2.4): 

dj - b + Mc^n, M > 1, (2.12) 

then all operators (2.3) will be precise as c — ► 0, whatever poor they may be for 
a big M and any fixed c > 0. On the contrary, the operator (2.3) with 

d. = d = max b. + c«y/n (2*13) 

J l<j<k J 

may be a “very good” one for any fixed c > 0 (e.g., if the difference max b. - 

min bj is small), but it is imprecise as c — ► 0. 

It is clear that the quality of “precision” describes the action of the 
distinction operator when X varies being always non-empty and the comparison 
cube becomes smaller as its edge c — ► 0. One cannot propose a distinction 
operator good for all cases. In each case an operator should be constructed 
according to the problem under consideration. 

3. THE BETA-ALGORITHM 

Starting with X, C as in (1.5), and N > 2, we denote B q = 0,1^ = {1, 2, 
..., N n } and proceed as follows. 

Iteration 1. Take x e X (not x e C as in /l/) and partition B = C 

0 v 0 ' 0 

into N n subcubes C }, such that c} n cf = </> for i ^ j and uc} = C. Basing 

on x q , as the representative of one of C?, say C.^, apply parallel translation of 

Cl.* to make it coincide, turn by turn, with each C.\ i = 1, 2, ..., N n ; then x 
10 i o 

e C } will define the representative x.^ e C.^ in each C. . This rule defines the 
10 1 1 i 

“translated grid generator” and the collection of x. , i = 1, 2, ..., N n , yields the 
grid for any particular choice of x q e C. 
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Compute all f(x?), i e 1^. 

Check the membership x } € X for each i e I q and define the index set 

J o = {i I x. 1 £ X, i £ y. (3.1) 

Clearly J q ^ <j> since x q e X. Compute 

P x = min f(x3) (3.2) 

ie J 

o 

and determine x^ such that f(x^) = p^. This constitutes the marginal 

comparison constant generator different from that employed in /l/. 

Delete all C? e B for which 
l o 

ffrj 1 ) - Pj > Fj, i £ I Q , Tj as in (1.7). (3.3) 

(This is a deletion operator analogous to that employed in /l/.) The remaining 

subcubes correspond to the index set: 

= 0 I ffrj 1 ) - Pj < r 1 , i £ y. (3.4) 

The closure of the subcubes with i e 1^ defines the set: 

Sj = {x I x £ c. 1 , i £ y, Sj ^ <j>. (3.5) 

Basing on x|, i c 1^, apply a precise distinction operator *■» $(X, C?) and 
exclude from further considerations every c} for which 4>. =» 0, i e L . The 

l l , 1 

remaining subcubes correspond to the index set: 

jj = {i | $. = i, i £ y, Jj ^ ^ (3.6) 

(note that x} e X cannot be deleted by this operation). The closure of the 

subcubes with i e defines the set: 
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Bj = {x I x € C.\ i £ Jj> c Sj c B o , Bj ^ 4 >- (3.7) 

Further iterations. Partition each C } e B. in the same way as B and 
1 1 o 

repeat Iteration 1 replacing x q , r^ by x? (ieJj), r^, etc., with r^ given by (1.7). 

In this process we come to the two monotonic sequences: 



p x > p 2 > ... > 


P > ... 
m 


(3.8) 


C = B D B, D B„ 3 
o — 1 — 2 ~ 


... D B D ... . 
“ m ” 


(3.9) 


Theorem 1. 


lim p = p = 
m 

m— >oo 


= min f(x) 
xeX 


(3.10) 


OO 


lim B = 0 B =? = 

m 1 ' m 


{x | f(x) = P°, X £ X}. 


(3.11) 



m— kx> m=0 



Proof . The proof is analogous to that given in /2/ with minor modifications 
due to the modifications in the algorithm. The principal points of the proof 
reflect the geometry of the action of the Beta-Algorithm and can be summarized 
as follows: 

(a) The construction of the deletion and distinction operators guarantees the 

O “O 

nonelimination of the global minimizers, that is, the points x e X . 

(b) The marginal comparison constant generator and successive partitions and 
deletions assure the monotonic descent within X . 

(c) By virtue of the Lemma for robust X and continuous f, the descent in (3.8) 

does not cease (although temporary stopovers may occur). This means that 

the process cannot get stuck at some p^ > p° because of the 
non-appearance of new basic points x^ yielding p^ = f(x^) < 
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(d) Due to (a)-(b)-(c), to the compacity of X and because -► 0 as m -*• oo, 

the limit p° is attained and also there exists lim x =» x° e X° c X. 

m 

m— >00 

(e) The precision property of the distinction operator guarantees that the entire 

o 

set X of all global minimizers and only of those minimizers is obtained in 
(3.11). □ 

Consider now the case when the distinction operator is imprecise or 
unavailable at all, the latter meaning that there is no exclusions by <£. = 0, so 

that J. = L in (3.6), = S t in (3.7) and S stand for B in (3.9) for m = 

ll vy l l v/ m m v/ 

1, 2, ... . 

Theorem 2 . 

If the Beta-Algorithm includes an imprecise distinction operator or no 
distinction operator, then 



lim p = p° = min f(x) (3.12) 

m — >00 m xeX 

lim x = x° 6 X° (3.13) 

m 

m— >-oo 



X 



oo 



Q„ B ' 



(3.14) 



Proof . The statement (3.12) is identical to (3.10) of Theorem 1 proved in 
(a)-(b)-(c)-(d) without the usage of the precision of the distinction operator, so it 
is correct. The statement (3.13) follows from the construction of the marginal 

comparison constant generator (3.2), from the ultimately non-ceasing descent in 
(3.8) and from the compactness of X. The statement (3.14) follows from the 
non-elimination of global minimizers. □ 



The Beta-Algorithm without a distinction operator is called Reduced 
Bet a- Algorithm. This algorithm is of slower convergence and does not generally 

Q 

determine the set X of all global minimizers but only some larger set containing 
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X°, see (3.14). 

4. RELATIONSHIP TO THE OPTIMALITY CONDITIONS 

Suppose that the compact robust set X in (1.1), (1.2) is given by the 
inequalities: 



X = {x 1 g.(x) > 0, i - 1, ..., m}. 

Consider the set X produced by the Bet a- Algorithm. If a point x 6 X is 
a global minimizer by virtue of certain sufficient conditions, then x° c X° due to 
nonelimination of global minimizers by the Beta- Algorithm. On the contrary, if 
x° is only a local minimizer, then x° / X° since all points which are not global 
minimizers are eliminated by the Bet a- Algorithm. Thus, relationship between the 
Bet a- Algorithm and sufficient optimality conditions is trivial. 

To obtain the necessary optimality conditions we have to impose the 
differentiability requirements on the functions involved and, in addition, certain 
constraint qualification. It can be demonstrated that currently used constraint 
qualifications imply the local robustness of X in the neighborhood of x° e X° c 
X, but clearly, not vice versa. This means that such requirements applied at a 
point x° e ? of the minimizing set X° produced by the Beta-Algorithm yield 
the necessary conditions corresponding to the qualifications imposed. The 
Bet a- Algorithm, however, does not need any constraint qualification except for the 
inf-robustness of X, that is, the robustness of X in the neighborhood of X . The 
algorithm is, thus, applicable to much larger class of problems than any particular 
method based on a set of necessary and/or sufficient conditions. It involves, 
however, an iterative process which may be avoided by application of necessary 
and sufficient conditions (at the cost of obtaining a single, maybe, local minimizer 
in nonconvex cases). 
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5. CONCLUSIONS 

We see that the Beta-Algorithm always finds the global minimum p° of f(x) 

over a compact robust set X and at least one of its global minimizers no matter 

whether it is equipped with a distinction operator of whatever quality. If it is 

equipped with a precise distinction operator, then it finds also the minimizing set 

X° = {x | f(x) = p°, x € X}. If the distinction operator is imprecise, then the 

oo 

algorithm delivers some larger set V — P| B containing X . 

m=0 

For the Cubic Algorithm /l/ distinction operators are redundant and the 
marginal comparison constant generator (3.2) coincides with the extremal generator 
proposed in /l/. 

The paper was written at the Chateau du Mont-Sainte-Anne, Beaupr6, 
Quebec, during the student ski-week. It is instructive that the idea of the 
semi-certain distinction operator and the principal lines of the algorithm came to 
the author in the chairlift on the way to the summit, Jan. 7-8, 1986. 
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ABSTRACT 

This paper contrasts the recent polynomial algorithms for linear programming 
of Khachian and Karmarkar. We show that each requires the solution of a 
weighted least-squares subproblem at every iteration. By comparing these 
subproblems we obtain further insights into the two methods. 

1. INTRODUCTION 

In this paper we contrast the ellipsoid method for linear inequalities, which 
was shown by Khachian [1979, 1980] to provide a polynomial algorithm for linear 
programming, with the new projective algorithm for linear programming of 
Karmarkar [1984a, b]. Note that the ellipsoid method was developed for convex, 
not necessarily differentiable, optimization by Yudin and Nemirovsky [1977] — see 
the survey article of Bland, Goldfarb and Todd [1981]. 

Many authors have noted that these two algorithms use many similar 
concepts — ideas from nonlinear programming, geometric motivation, and infinite 
iterative schemes that can be truncated after a polynomial number of steps when 
applied to rational data, with an exact optimal solution then available by 
rounding. However, the details of the two algorithms seem very different. We 
will show that in fact the heart of each iteration of either algorithm is the 
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solution of a weighted least-squares subproblem, and that these subproblems are 
very closely related. This viewpoint allows further insights into the two methods, 
in particular suggesting reasons for the very slow convergence of the ellipsoid 
method compared to the apparently very fast convergence of the projective 
algorithm. 

The weighted least-squares subproblems have other important features. Both 
the ellipsoid and the projective algorithms appear at first sight not to provide 
solutions to the dual linear programming problem, but a closer examination shows 
that dual solutions are indeed generated during the course of the methods, 
essentially from the least-squares problems. Naturally, optimal dual solutions can 
be generated from optimal primal solutions at termination; however, approximate 
dual solutions at intermediate stages are very useful in guaranteeing the quality of 
the current solution or in certifying infeasiblity. 

Section 2 describes the ellipsoid method from the viewpoint of Burrell and 
Todd [1985], showing how the weighted least-squares subproblems arise. Section 3 
then outlines Karmarkar’s algorithm, demonstrating that computation of the search 
direction once again requires the solution of a least-squares subproblem. In 
section 4 we contrast the two methods by comparing these subproblems and how 
they change from one iteration to the next. 

In the rest of this introductory section we discuss the two methods 
informally. As we have noted above, for real data both algorithms are infinite 
iterative methods, for which it is possible to establish linear convergence rates. 
Such behavior is usually regarded as unacceptably slow in nonlinear programming, 
where under suitable conditions (quasi-) Newton methods achieve superlinear rates; 
however, this linear convergence rate is from the very first iteration, rather than 
asymptotic, and the convergence ratio can be bounded as a function of the 
dimensions alone, independent of the data. 

Using this ratio enables one to establish polynomial bounds on the number of 
bit operations required to solve a linear programming problem. Suppose a 

problem has n variables, a number of constraints of the same order of magnitude, 
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and input size L (i.e., L is the number of bits to specify the data of the 

4 2. 

problem). Then the ellipsoid method requires about 0(n L ) bit operations and 

3 5 2 

Karmarkar’s algorithm about 0(n * L ). Note that Karmarkar [1984a, b] quotes 
ft 2 

0(n L ) for the ellipsoid method, from Khachian [1979], whereas the revised figure 

appears in Khachian [1980]. Also, the basic algorithm in Karmarkar [1984a, b] 
4 2 

requires about 0(n L ) operations — this is the version we shall describe in 

5 

section 3 — whereas the n* factor is removed by a modification that reduces the 
linear algebra necessary at each step. Thus there is only a slight advantage to 
Karmarkar’s algorithm from the standpoint of theoretical bounds. 

In practical performance, the differences become marked. The ellipsoid 
method appears to require a number of iterations close to its worst-case bound 
0(n L). On the other hand, a number of studies have established that many 
variants of Karmarkar’s algorithm only take a number of iterations between 20 
and 50 to get a very accurate solution, and this number appears to grow very 

slowly, if at ail, with n. It is not our aim to discuss the computational 

significance of Karmarkar’s algorithm, but at least for certain very large sparse 
problems exceptional results have been obtained. Besides the conference 

presentations of Karmarkar and Karmarkar and Sinha [1985], the best 

computational results have been given by Gill, Murray, Saunders, Tomlin and 
Wright [1985] and Adler, Karmarkar, Resende and Veiga [1986]. The latter paper 
reports times ranging from three times slower to eight times faster than the 

MINOS code on a set of medium to large problems. 

In contrast to the overwhelming computational superiority of the projective 
algorithm over the ellipsoid method, its theoretical implications are far more 
limited. Grotschel, Lov&sz and Schrijver [1981, 1986] and others have used the 
ellipsoid method to show that certain combinatorial optimization problems are in 
P (i.e., have polynomial algorithms) and others are NP-hard (i.e., are unlikely to 
have polynomial algorithms). A key to its use is the fact that the ellipsoid 
method does not need to have all the constraints listed in advance — they can 
be generated as needed. Thus it is possible to handle problems with an 

exponential number of constraints. On the other hand, it appears that 
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Karmarkar’s algorithm requires all the constraints and variables to be explicitly 
present, and thus it cannot be used to solve such problems with our present 
knowledge. 

2. THE ELLIPSOID METHOD 

The preferred problem for the ellipsoid method is that of determining 
feasibility of a system of linear inequalities. To facilitate comparison with 
Karmarkar’s algorithm, we assume that we seek a point in 

Y - {y e IR m : A T y < c} (2.1) 

where A is an mxn matrix and c an n- vector. We assume that a bound is 
known for all solutions to Y: 

(Al) There is a known y^ e lR m and R > 0 such that Y c S(y^, R), 

where S(y°, R) denotes the ball {y e IR m : ||y - y°|| < R} and ||-|| denotes the 
Euclidean norm. 

The ellipsoid method generates a sequence of ellipsoids {E^} with centers 
{yk} such that 

E q = S(y°, R); (2.2) 

Y c E k , k > 0; (2.3) 

if y k ^ Y, then 

Yo1 E k+l/ Vo1 E k < ex P( _1 / 2 ( m + !)) • ( 2 - 4 ) 

Of course, if y e Y for some k, the algorithm stops; otherwise by (2.4) the 
ellipsoids generated have volumes shrinking geometrically to zero. If we further 



assume 
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(A2) If Y ^ there is some y e !R m , p > 0, such that S(£, p) c Y, 

k 

then if the algorithm does not generate a feasible iterate y in 

[2m(m + l)tn(R/p)] + 1 steps, we can stop and conclude that Y = <t>. This 

follows from (2.2)-(2.4); in this number of steps, the volume of has shrunk 
from that of a ball of radius R to less than that of one of radius p, and by 
(A2) Y must be empty. Assumptions (Al) and (A2) imply that Y is bounded 
and has a nonempty interior; if (Al) is not satisfied, the algorithm cannot be 
started while if (A2) fails it will usually not terminate. If the data in A and c 
are integer with an input size L, then bounds can be added to the variables and 
the right-hand sides c increased a little so that the resulting polyhedron Y 
satisfies (Al) and (A2). Moreover, Y will be nonempty iff Y is for suitable 
bounds (of order 2^) and perturbations of c (of order 2 ^). In this way, one 
can show that solving the feasibility problem for a polytope Y satisfying (Al) and 
(A2) in polynomial time allows a polynomial algorithm for general linear 
programming. For more details, see for example Bland, Goldfarb and Todd 
[1981]. 

k 

The ellipsoid is usually represented by its center y and a symmetric 
positive definite matrix 

E fc = {y e R m : (y - y k ) T B fc 1 (y - y k ) < 1} • (2-5) 

Ignoring finite precision, we can give the update formulae for E^^ ^ foll° ws - If 

i T k 

y £ Y, find some violated constraint a y < 7 in (2.1), so that ay > 7- 

Then 

y k+l = y k _ r k B k a/(a T B t a) a , (2-6) 

B _ , ( B , 

B k + i - M k - "k t J - 



(2.7) 
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The simplest method to choose the scalar parameters is 

r k = !/(m + 1). * k = 2/(m + 1) , = m 2 /(m 2 - 1), 

T k 

but different values depending on a y -7 can give a greater decrease in the 

volume of E, , . 

k+1 

One of the disadvantages of this implementation of the ellipsoid method is 

that it provides no easy way to prove infeasibility. Usually, if Y =* </>, at a 

certain k we will have 

{y e E fc : a T y < 7} = <j> , (2.8) 

T 

where a y < 7 is one of the inequalities in (2.1). However, if round-off errors 

may have occurred during the algorithm, we cannot be certain that (2.3) holds. 

A much more satisfactory way to demonstrate infeasibility would be to generate a 
vector x e H n with 

T 

Ax = 0, c x < 0, x > 0 . (2.9) 

i.e., a solution to the alternative system in Parkas’ lemma. Similarly, if we are 
using the ellipsoid method to solve 

0 T 

min{(a ) y : y e Y} . (2.10) 

then we would like to use the value 

min{(a°) T y : y 6 E k > (2.11) 

to provide a lower bound guaranteeing the quality of a feasible solution. 

However, in the presence of round-off error, such a bound may not be valid. 

Linear programming duality provides a more satisfactory way to derive a bound: 
T 

— c x is a lower bound if x e 1R satisfies 



Ax == -a^ , x > 0. 



( 2 . 12 ) 




